GVS | B Opleiding tot Verpleegkundige (HBO-V) - voltijd - versie 1.0
Auteur
Theo Bakker, lector Learning Technology & Analytics, De HHs
Publicatiedatum
6 juni 2024
1 Inleiding
1.1 Het nut van prognosemodellen
Prognosemodellen kunnen inzicht bieden in de factoren die gecorreleerd zijn aan de uitval van studenten. Met deze inzichten kan een opleiding interventies ontwikkelen om de uitval te verminderen of te voorkomen. Voor het lectoraat is het van belang dat we in vervolg hierop kunnen analyseren of er sprake is van bias in de data in relatie tot uitval en er mogelijk sprake is van een gebrek aan kansengelijkheid. Een prognosemodel is dus de opmaat naar een fairness analyse. Zie voor een verdere toelichting het onderzoeksprogramma ‘No Fairness without Awareness’.
1.2 Toelichting op de methode
Voor de ontwikkeling van prognosemodellen gebruiken we de aanpak van Tidymodels. Tidymodels is een framework voor het bouwen van een prognosemodel. Hiermee verzekeren we ons van een systematische, herhaalbare en schaalbare aanpak voor het bouwen van prognosemodellen.
1.3 Toelichting op de data
De basis voor deze analyse is studiedata van De Haagse Hogeschool (De HHs), verrijkt door het lectoraat Learning Technology & Analytics. De data bevat informatie over de inschrijvingen van studenten in het eerste jaar van de opleiding:
Demografische kenmerken: geslacht, leeftijd, reistijd en SES totaalscore.
Vooropleidingskenmerken: toelaatgevende vooropleiding, studiekeuzeprofiel, gemiddeld eindcijfer in de vooropleiding en eventuele deelname aan het Navitas programma.
Aanmeldingskenmerken: aansluiting (direct na diploma, tussenjaar, switch), dag van aanmelding, aantal parallelle studies aan De HHs en collegejaar.
1.4 Toelichting op de analyse
We toetsen in deze analyse Uitval na 1 jaar bij studenten zonder propedeusediploma, voortaan Uitval genoemd.
Uitval is gedefinieerd als het niet meer ingeschreven staan in de opleiding in een aansluitend collegejaar. Een wisseling van opleidingsvorm binnen de opleiding, bijv. van voltijd in jaar 1 naar duaal in jaar 2, geldt niet als uitval.
2 Voorbereidingen
2.1 Laad de data
We laden een subset in van historische data specifiek voor:
Opleiding: GVS | B Opleiding tot Verpleegkundige (HBO-V), voltijd, eerstejaars - Uitval na 1 jaar
Toon code
## Laad de data voor de HBO-V opleidingdfOpleiding_inschrijvingen_base <-get_lta_studyprogram_enrollments_pin(board ="HHs/Inschrijvingen",faculty = faculteit,studyprogram = opleidingsnaam_huidig,studytrack = opleiding,studyform =toupper(opleidingsvorm),range ="eerstejaars")## Herschik de levelsSet_Levels(dfOpleiding_inschrijvingen_base)dfOpleiding_inschrijvingen_base <- dfOpleiding_inschrijvingen_base |>## Maak een eenvoudige Uitval variabele aanMutate_Uitval(sUitval_model) |>## Maak van de uitval variabele een factormutate(SUC_Uitval =as.factor(SUC_Uitval)) |>## Verbijzonder eventueel op basis van het propedeusediplomaFilter_Propedeusediploma(sPropedeusediploma) |>## Maak van de Dubbele studie variabele een Ja/Nee variabelemutate(INS_Dubbele_studie =ifelse(INS_Aantal_inschrijvingen >1, "Ja", "Nee")) |>## Verwijder INS_Aantal_inschrijvingenselect(-INS_Aantal_inschrijvingen) |>## Pas voor een aantal variabelen de levels aanMutate_Levels(c("VOP_Studiekeuzeprofiel_LTA_afkorting","INS_Aansluiting_LTA","VOP_Toelaatgevende_vooropleiding_soort" ),list(lLevels_skp, lLevels_vop, lLevels_vop) )## B Huidtherapie: Filter op uitsluitend studenten met een rangnummer (selectie)if(opleiding =="HDT") { dfOpleiding_inschrijvingen_base <- dfOpleiding_inschrijvingen_base |>filter(!is.na(RNK_Rangnummer)) }
2.2 Selecteer en inspecteer de data
We selecteren eerst de relevante variabelen. We verwijderen daarbij variabelen die maar 1 waarde hebben. We bekijken de variabelen in een samenvatting in relatie tot de Uitval. Daarnaast bekijken we de kwaliteit van de data op missende waarden.
Toon code
lSelect <-c("INS_Student_UUID_opleiding_vorm","CBS_APCG_tf","DEM_Geslacht","DEM_Leeftijd_1_oktober","GIS_Tijd_fiets_OV","INS_Collegejaar","INS_Dagen_tussen_aanmelding_en_1_september","INS_Dubbele_studie","INS_Aansluiting_LTA","INS_Navitas_tf","SES_Deelscore_arbeid","SES_Deelscore_welvaart","SES_Totaalscore","SUC_Uitval","VOP_Gemiddeld_cijfer_cijferlijst","VOP_Gemiddeld_eindcijfer_VO_van_de_hoogste_vooropleiding_voor_het_HO","VOP_Cijfer_CE1_nederlands","VOP_Cijfer_CE1_engels","VOP_Cijfer_CE_proxy_wiskunde","VOP_Cijfer_CE1_natuurkunde","VOP_Studiekeuzeprofiel_LTA_afkorting","VOP_Toelaatgevende_vooropleiding_soort" )## B Huidtherapie: voeg de variabele RNK_Rangnummer toeif(opleiding =="HDT") { lSelect <-c(lSelect, "RNK_Rangnummer")}## Maak een subsetdfOpleiding_inschrijvingen <- dfOpleiding_inschrijvingen_base |>## Selecteer de relevante variabelenselect_at(lSelect) |>## Hernoem variabelen voor beter leesbare namenrename(ID = INS_Student_UUID_opleiding_vorm,Geslacht = DEM_Geslacht,Leeftijd = DEM_Leeftijd_1_oktober,Reistijd = GIS_Tijd_fiets_OV,Dubbele_studie = INS_Dubbele_studie,Collegejaar = INS_Collegejaar,Aanmelding = INS_Dagen_tussen_aanmelding_en_1_september,Aansluiting = INS_Aansluiting_LTA,Navitas = INS_Navitas_tf,APCG = CBS_APCG_tf,SES_Arbeid = SES_Deelscore_arbeid,SES_Welvaart = SES_Deelscore_welvaart,SES_Totaal = SES_Totaalscore, Uitval = SUC_Uitval,Cijfer_SE_VO = VOP_Gemiddeld_cijfer_cijferlijst,Cijfer_CE_VO = VOP_Gemiddeld_eindcijfer_VO_van_de_hoogste_vooropleiding_voor_het_HO,Cijfer_CE_Nederlands = VOP_Cijfer_CE1_nederlands,Cijfer_CE_Engels = VOP_Cijfer_CE1_engels,Cijfer_CE_Wiskunde = VOP_Cijfer_CE_proxy_wiskunde,Cijfer_CE_Natuurkunde = VOP_Cijfer_CE1_natuurkunde,Studiekeuzeprofiel = VOP_Studiekeuzeprofiel_LTA_afkorting,Vooropleiding = VOP_Toelaatgevende_vooropleiding_soort ) |>## Pas CBS_APCG_tf aan naar factormutate(APCG =case_when(APCG ==TRUE~"Ja", APCG ==FALSE~"Nee",.default ="Onbekend")) |>## Geef aan waar missende cijfers in het VO zijnMutate_Cijfers_VO() |>## Verwijder variabelen, waarbij er maar 1 waarde isselect(where(~n_distinct(.) >1)) |>arrange(Collegejaar, ID)## B Huidtherapie: hernoem de variabele RNK_Rangnummerif(opleiding =="HDT") { dfOpleiding_inschrijvingen <- dfOpleiding_inschrijvingen |>rename(Rangnummer = RNK_Rangnummer)} dfOpleiding_inschrijvingen <- dfOpleiding_inschrijvingen |> ltabase::sort_distinct()## Verwijder de basis datasetrm(dfOpleiding_inschrijvingen_base)
Studentkenmerken versus Uitval
Variabele
Uitval
p-value2
Totaal, N = 18181
Ja, N=727 (40%)1
Nee, N=1091 (60%)1
Aanmelding
133,20 (74,18)
153,98 (72,88)
<0,001
145,67 (74,08)
Aansluiting
0,025
2e Studie
6 (32%)
13 (68%)
19 (100%)
Direct
357 (38%)
579 (62%)
936 (100%)
Na CD
4 (16%)
21 (84%)
25 (100%)
Overig
2 (67%)
1 (33%)
3 (100%)
Switch extern
230 (42%)
314 (58%)
544 (100%)
Switch intern
37 (51%)
36 (49%)
73 (100%)
Tussenjaar
91 (42%)
127 (58%)
218 (100%)
APCG
<0,001
Ja
305 (46%)
364 (54%)
669 (100%)
Nee
342 (35%)
629 (65%)
971 (100%)
Onbekend
80 (45%)
98 (55%)
178 (100%)
Cijfer_CE_Engels
6,84 (1,13)
6,47 (1,22)
<0,001
6,60 (1,21)
Cijfer_CE_Engels_missing
<0,001
Ja
438 (45%)
539 (55%)
977 (100%)
Nee
289 (34%)
552 (66%)
841 (100%)
Cijfer_CE_Natuurkunde
6,04 (1,00)
6,21 (0,97)
0,23
6,15 (0,98)
Cijfer_CE_Natuurkunde_missing
0,016
Ja
637 (41%)
911 (59%)
1.548 (100%)
Nee
90 (33%)
180 (67%)
270 (100%)
Cijfer_CE_Nederlands
6,12 (0,92)
6,17 (0,85)
0,45
6,15 (0,88)
Cijfer_CE_Nederlands_missing
<0,001
Ja
437 (45%)
539 (55%)
976 (100%)
Nee
290 (34%)
552 (66%)
842 (100%)
Cijfer_CE_VO
6,38 (0,37)
6,52 (0,43)
<0,001
6,48 (0,41)
Cijfer_CE_VO_missing
<0,001
Ja
414 (48%)
455 (52%)
869 (100%)
Nee
313 (33%)
636 (67%)
949 (100%)
Cijfer_CE_Wiskunde
6,47 (1,10)
6,49 (1,18)
0,95
6,48 (1,16)
Cijfer_CE_Wiskunde_missing
<0,001
Ja
466 (45%)
578 (55%)
1.044 (100%)
Nee
261 (34%)
513 (66%)
774 (100%)
Cijfer_SE_VO
6,39 (0,36)
6,52 (0,43)
<0,001
6,48 (0,41)
Cijfer_SE_VO_missing
<0,001
Ja
430 (46%)
502 (54%)
932 (100%)
Nee
297 (34%)
589 (66%)
886 (100%)
Dubbele_studie
0,066
Ja
13 (59%)
9 (41%)
22 (100%)
Nee
714 (40%)
1.082 (60%)
1.796 (100%)
Geslacht
<0,001
M
177 (56%)
140 (44%)
317 (100%)
V
550 (37%)
951 (63%)
1.501 (100%)
Leeftijd
20,75 (3,43)
19,95 (3,08)
<0,001
20,27 (3,25)
Reistijd
32,18 (18,05)
32,07 (16,62)
0,50
32,11 (17,20)
SES_Arbeid
-0,03 (0,10)
-0,01 (0,10)
<0,001
-0,02 (0,10)
SES_Totaal
-0,14 (0,35)
-0,07 (0,34)
<0,001
-0,10 (0,35)
SES_Welvaart
-0,06 (0,16)
-0,03 (0,16)
<0,001
-0,04 (0,16)
Studiekeuzeprofiel
<0,001
EM
51 (40%)
78 (60%)
129 (100%)
CM
54 (38%)
90 (63%)
144 (100%)
EM&CM
27 (25%)
81 (75%)
108 (100%)
NT
12 (80%)
3 (20%)
15 (100%)
NG
175 (35%)
318 (65%)
493 (100%)
NT&NG
50 (29%)
120 (71%)
170 (100%)
OS
0 (0%)
3 (100%)
3 (100%)
CERT
0 (0%)
2 (100%)
2 (100%)
ALG
3 (50%)
3 (50%)
6 (100%)
BI
1 (100%)
0 (0%)
1 (100%)
EA
61 (61%)
39 (39%)
100 (100%)
HO
9 (38%)
15 (63%)
24 (100%)
HB
6 (55%)
5 (45%)
11 (100%)
ICT
0 (0%)
3 (100%)
3 (100%)
MedV
8 (47%)
9 (53%)
17 (100%)
TP
1 (20%)
4 (80%)
5 (100%)
TR
6 (40%)
9 (60%)
15 (100%)
TSL
2 (67%)
1 (33%)
3 (100%)
UV
6 (50%)
6 (50%)
12 (100%)
VS
2 (33%)
4 (67%)
6 (100%)
VNL
6 (60%)
4 (40%)
10 (100%)
ZW
173 (52%)
160 (48%)
333 (100%)
Vooropleiding
<0,001
MBO
286 (52%)
262 (48%)
548 (100%)
HAVO
344 (35%)
636 (65%)
980 (100%)
VWO
25 (30%)
59 (70%)
84 (100%)
BD
47 (46%)
56 (54%)
103 (100%)
HO
19 (26%)
53 (74%)
72 (100%)
CD
6 (19%)
25 (81%)
31 (100%)
1 Mean (SD); n (%)
2 Wilcoxon rank sum test; Fisher’s Exact Test for Count Data with simulated p-value (based on 2000 replicates); Pearson’s Chi-squared test
Toon code
## Toon een samenvatting van de data, gesorteerd op missende waardendiagnose(dfOpleiding_inschrijvingen) |>mutate(missing_percent =round(missing_percent, 2),unique_rate =round(missing_percent, 2)) |>arrange(desc(missing_percent)) |> knitr::kable(caption ="Kwaliteit van de data voor bewerkingen (gesorteerd op missende waarden)")
Kwaliteit van de data voor bewerkingen (gesorteerd op missende waarden)
variables
types
missing_count
missing_percent
unique_count
unique_rate
Cijfer_CE_Natuurkunde
numeric
1548
85.15
47
85.15
Cijfer_CE_Wiskunde
numeric
1044
57.43
61
57.43
Cijfer_CE_Engels
numeric
977
53.74
66
53.74
Cijfer_CE_Nederlands
numeric
976
53.69
52
53.69
Cijfer_SE_VO
numeric
932
51.27
27
51.27
Cijfer_CE_VO
numeric
869
47.80
27
47.80
Studiekeuzeprofiel
factor
208
11.44
23
11.44
SES_Arbeid
numeric
178
9.79
311
9.79
SES_Totaal
numeric
178
9.79
592
9.79
SES_Welvaart
numeric
178
9.79
421
9.79
Reistijd
numeric
29
1.60
379
1.60
Aanmelding
numeric
0
0.00
319
0.00
Aansluiting
factor
0
0.00
7
0.00
APCG
character
0
0.00
3
0.00
Cijfer_CE_Engels_missing
character
0
0.00
2
0.00
Cijfer_CE_Natuurkunde_missing
character
0
0.00
2
0.00
Cijfer_CE_Nederlands_missing
character
0
0.00
2
0.00
Cijfer_CE_VO_missing
character
0
0.00
2
0.00
Cijfer_CE_Wiskunde_missing
character
0
0.00
2
0.00
Cijfer_SE_VO_missing
character
0
0.00
2
0.00
Collegejaar
numeric
0
0.00
11
0.00
Dubbele_studie
character
0
0.00
2
0.00
Geslacht
character
0
0.00
2
0.00
ID
character
0
0.00
1818
0.00
Leeftijd
integer
0
0.00
23
0.00
Uitval
factor
0
0.00
2
0.00
Vooropleiding
factor
0
0.00
6
0.00
2.3 Bewerk de data
Uit de eerste diagnose blijkt dat niet alle variabelen goed genoeg zijn voor het bouwen van een prognosemodel: er zijn missende waarden en niet alle veldtypes zijn geschikt. We passen de variabelen aan zodat we in het model er goed mee kunnen werken.
Prognosemodellen kunnen niet omgaan met missende waarden. Om bias te voorkomen verwijderen we geen rijen met missende waarden, maar vullen die op (imputatie). We bewerken de data zo dat alle missende waarden worden opgevuld: bij numerieke waarden met het gemiddelde en bij categorische variabelen met ‘Onbekend’.
We passen sommige variabelen aan, zodat ze in het model gebruikt kunnen worden: tekstvelden zetten we om naar factor; logische variabelen (Ja/Nee) zetten we om naar een numerieke variabele (1/0).
De uitkomstvariabele, Uitval, leiden we af van de variabele SUC_Uitval_aantal_jaar_LTA. Als de waarde daar 1 is, is de student na 1 jaar uitgevallen, 2 na 2 jaar, etc.
Een fictief studentnummer (INS_Student_UUID_opleiding_vorm) gebruiken we ook, zodat we - als er afwijkende resultaten zijn - de dataset gericht kunnen onderzoeken indien nodig.
Toon code
## Bewerk de datadfOpleiding_inschrijvingen <- dfOpleiding_inschrijvingen |>## Imputeer alle numerieke variabelen met de meanmutate(across(where(is.numeric), ~ifelse(is.na(.x),mean(.x, na.rm = T), .x )) ) |>## Zet character variabelen om naar factormutate(across(where(is.character), as.factor)) |>## Zet logische variabelen om naar 0 of 1mutate(across(where(is.logical), as.integer)) |>## Vul in factoren missende waarden op met "Onbekend"mutate(across(where(is.factor), ~suppressWarnings(fct_explicit_na(.x, na_level ="Onbekend") ))) |>## Herschik de kolommen, zodat Uitval vooraan staatselect(Uitval, everything()) ## Bekijk de data## glimpse(dfOpleiding_inschrijvingen) ## Maak een diagnose van de datadiagnose(dfOpleiding_inschrijvingen) |>mutate(missing_percent =round(missing_percent, 2),unique_rate =round(unique_rate, 2)) |> knitr::kable(caption ="Kwaliteit van de data na bewerkingen")
Kwaliteit van de data na bewerkingen
variables
types
missing_count
missing_percent
unique_count
unique_rate
Uitval
factor
0
0
2
0.00
Aanmelding
numeric
0
0
319
0.18
Aansluiting
factor
0
0
7
0.00
APCG
factor
0
0
3
0.00
Cijfer_CE_Engels
numeric
0
0
66
0.04
Cijfer_CE_Engels_missing
factor
0
0
2
0.00
Cijfer_CE_Natuurkunde
numeric
0
0
47
0.03
Cijfer_CE_Natuurkunde_missing
factor
0
0
2
0.00
Cijfer_CE_Nederlands
numeric
0
0
52
0.03
Cijfer_CE_Nederlands_missing
factor
0
0
2
0.00
Cijfer_CE_VO
numeric
0
0
27
0.01
Cijfer_CE_VO_missing
factor
0
0
2
0.00
Cijfer_CE_Wiskunde
numeric
0
0
61
0.03
Cijfer_CE_Wiskunde_missing
factor
0
0
2
0.00
Cijfer_SE_VO
numeric
0
0
27
0.01
Cijfer_SE_VO_missing
factor
0
0
2
0.00
Collegejaar
numeric
0
0
11
0.01
Dubbele_studie
factor
0
0
2
0.00
Geslacht
factor
0
0
2
0.00
ID
factor
0
0
1818
1.00
Leeftijd
integer
0
0
23
0.01
Reistijd
numeric
0
0
379
0.21
SES_Arbeid
numeric
0
0
311
0.17
SES_Totaal
numeric
0
0
592
0.33
SES_Welvaart
numeric
0
0
421
0.23
Studiekeuzeprofiel
factor
0
0
23
0.01
Vooropleiding
factor
0
0
6
0.00
2.4 Bekijk de onderlinge correlaties
Het is verstandig om voorafgaand aan het bouwen van een model te kijken naar de onderlinge correlaties tussen numerieke variabelen. Dit geeft inzicht in de data en kan helpen bij het maken van keuzes voor het model of de duiding van de uitkomsten.
Toon code
## Maak een plot van de onderlinge correlaties in numerieke variabelendfOpleiding_inschrijvingen |>select(-Collegejaar) |>select(where(is.numeric)) |>cor() |> corrplot::corrplot(order ='hclust', addrect =4,method ="number", tl.cex =0.8, tl.col ="black",diag =FALSE)
2.5 Bouw de trainingset, testset en validatieset
De data is nu geschikt om een prognosemodel mee te bouwen.
Om het model te bouwen, testen en valideren, splitsen we de data in drie delen van 60%, 20% en 20%. We doen dit op zo’n manier, dat elk deel ongeveer een gelijk aantal studenten bevat dat uitvalt.
We trainen het model op basis van 60% en testen de modellen tijden het trainen op de overige 20% (de testset).
Als het model klaar is, valideren we het op de 20% studenten uit de validatieset, die is opgebouwd uit zo recent mogelijke data. De validatieset blijft dus de gehele tijd ongemoeid, zodat we overfitting - een te goed model op bekende data, maar slechte performance op onbekende data - voorkomen.
Een willekeurig, maar vaststaand seedgetal voorkomt dat we bij elke run van het model c.q. deze code een net iets andere uitkomst krijgen.
Toon code
set.seed(0821)## Splits de data in 3 delen: 60%, 20% en 20%splits <-initial_validation_split(dfOpleiding_inschrijvingen,strata = Uitval,prop =c(0.6, 0.2))## Maak drie sets: een trainingset, een testset en een validatiesetdfUitval_train <-training(splits)dfUitval_test <-testing(splits)dfUitval_validation <-validation_set(splits)## Maak een resample set op basis van 10 folds (default)dfUitval_resamples <-vfold_cv(dfUitval_train, strata = Uitval)## Training set proporties uitvaldfUitval_train |>count(Uitval) |>mutate(prop = n/sum(n)) |> knitr::kable(caption ="Trainingsset")
Trainingsset
Uitval
n
prop
FALSE
654
0.6
TRUE
436
0.4
Toon code
## Test set proporties uitvaldfUitval_test |>count(Uitval) |>mutate(prop = n/sum(n)) |> knitr::kable(caption ="Testset")
Testset
Uitval
n
prop
FALSE
219
0.6
TRUE
146
0.4
3 Model I: Logistische Regressie
Het eerste model is een logistische regressie met penalized likelihood; we gebruiken de glmnet engine voor het bouwen van het model. Penalized likelihood is een techniek die helpt bij het voorkomen van overfitting. Glmnet is een populair package voor het bouwen van logistische regressiemodellen.
We gebruiken de Area under the ROC Curve (AUC/ROC) als performance metric.
Vervolgens zetten we meerdere stappen in een ‘recipe’. We verwijderen de student ID en het collegejaar uit de data, omdat deze niet moet worden gebruikt in het model. We converteren factoren naar dummy variabelen, verwijderen zero values en centreren en schalen numerieke variabelen.
## Bouw de recipe: logistische regressielr_recipe <-recipe(Uitval ~ ., data = dfUitval_train) |>update_role(ID, new_role ="ID") |>## Zet de student ID als ID variabelestep_rm(ID, Collegejaar) |>## Verwijder ID en collegejaar uit het modelstep_dummy(all_nominal_predictors()) |>## Converteer factoren naar dummy variabelenstep_zv(all_predictors()) |>## Verwijder zero valuesstep_normalize(all_numeric_predictors()) ## Centreer en schaal numerieke variabelen## Toon de recipetidy(lr_recipe) |> knitr::kable()
number
operation
type
trained
skip
id
1
step
rm
FALSE
FALSE
rm_yTMrb
2
step
dummy
FALSE
FALSE
dummy_0LqEE
3
step
zv
FALSE
FALSE
zv_P3EsP
4
step
normalize
FALSE
FALSE
normalize_Yy34e
3.3 Maak de workflow
Voor de uitvoering bouwen we een nieuwe workflow. Daaraan voegen we het model en de bewerkingen in de recipe toe.
## Maak de workflow: logistische regressielr_workflow <-workflow() |>## Maak een workflowadd_model(lr_mod) |>## Voeg het model toeadd_recipe(lr_recipe) ## Voeg de recipe toe## Toon de workflowlr_workflow
Het model moet getuned worden. Dit houdt in dat we de beste parameters voor het model moeten vinden. We maken een grid met verschillende penalty waarden. Daarmee kunnen we vervolgens het beste model selecteren met de hoogste ROC/AUC. We plotten de resultaten van de tuning, zodat we hieruit het beste model kunnen kiezen.
## Maak een grid: logistische regressielr_reg_grid <-tibble(penalty =10^seq(-4, -1, length.out =30))## Train en tune het model: logistische regressielr_res <- lr_workflow |>tune_grid(dfUitval_validation,grid = lr_reg_grid,control =control_grid(save_pred =TRUE),metrics =metric_set(roc_auc))
Toon code
## Plot de resultaten + een rode verticale lijn voor de max AUClr_plot <- lr_res |>collect_metrics() |>ggplot(aes(x = penalty, y = mean)) +geom_point() +geom_line() +ylab("Area under the ROC Curve") +scale_x_log10(labels = scales::label_number())# Zoek de penalty waarde met de max AUCmax_auc_penalty <- lr_res |>collect_metrics() |>filter(mean ==max(mean)) |>pull(penalty)# Voeg de rode verticale lijn toe aan lr_plotlr_plot_plus <- lr_plot +geom_vline(xintercept = max_auc_penalty, color ="red")# Vindt een mean voor de max AUC die hoger ismax_auc_mean <- lr_res |>collect_metrics() |>filter(mean ==max(mean)) |>pull(penalty)lr_plot_plus
3.5 Kies het beste model
We evalueren modellen met een zo hoog mogelijke Area under the ROC Curve (AUC/ROC) en een zo laag mogelijke penalty. Zo kunnen we uit de resultaten het beste model kiezen. Tot slot maken we een ROC curve om de prestaties van het model te visualiseren.
## Toon het beste modeltop_models <- lr_res |>show_best(metric ="roc_auc", n =10) |>mutate(mean =round(mean, 6)) |>arrange(penalty) top_models|> knitr::kable()
penalty
.metric
.estimator
mean
n
std_err
.config
0.0035622
roc_auc
binary
0.652199
1
NA
Preprocessor1_Model16
0.0045204
roc_auc
binary
0.652958
1
NA
Preprocessor1_Model17
0.0057362
roc_auc
binary
0.654160
1
NA
Preprocessor1_Model18
0.0072790
roc_auc
binary
0.654919
1
NA
Preprocessor1_Model19
0.0092367
roc_auc
binary
0.655489
1
NA
Preprocessor1_Model20
0.0117210
roc_auc
binary
0.655394
1
NA
Preprocessor1_Model21
0.0148735
roc_auc
binary
0.653243
1
NA
Preprocessor1_Model22
0.0188739
roc_auc
binary
0.655220
1
NA
Preprocessor1_Model23
0.0239503
roc_auc
binary
0.656675
1
NA
Preprocessor1_Model24
0.0303920
roc_auc
binary
0.655315
1
NA
Preprocessor1_Model25
## Selecteer het beste model: logistische regressielr_best <- lr_res |>collect_metrics() |>filter(mean ==max(mean)) |>slice(1) lr_best|>mutate(mean =round(mean, 6)) |> knitr::kable()
penalty
.metric
.estimator
mean
n
std_err
.config
0.0239503
roc_auc
binary
0.656675
1
NA
Preprocessor1_Model24
## Verzamel de predicties en evalueer het model (AUC/ROC): logistische regressielr_auc <- lr_res |>collect_predictions(parameters = lr_best) |>roc_curve(Uitval, .pred_FALSE) |>mutate(model ="Logistisch Regressie")autoplot(lr_auc)
## Bepaal de AUC van het beste modellr_auc_highest <- lr_res |>collect_predictions(parameters = lr_best) |>roc_auc(Uitval, .pred_FALSE)## Voeg de naam van het model en de AUC toe dfModel_resultsdfModel_results <- dfModel_results |>add_row(model ="Logistic Regression", auc = lr_auc_highest$.estimate)
4 Model II: Tree-based ensemble
Het tweede model is een random forest: een ensemble van beslisbomen (decision trees). Het is een krachtig model dat goed om kan gaan met complexe data en veel variabelen.
We gebruiken de ranger engine voor het bouwen van het model.
4.1 Bepaal het aantal PC-cores
Omdat een random forest model veel berekeningen vereist, willen we daarvoor alle computerkracht gebruiken die beschikbaar is. Het aantal CPU’s (cores) van de computer bepaalt hoe snel het model getraind kan worden. Deze informatie gebruiken we bij het bouwen van het model.
Toon code
## Bepaal het aantal corescores <- parallel::detectCores()
4.2 Maak het model
We bouwen eerst het model. We gebruiken de rand_forest functie om het model te bouwen. We tunen de mtry en min_n parameters. De mtry parameter bepaalt het aantal variabelen dat per boom wordt gebruikt. De min_n parameter bepaalt het minimum aantal observaties dat in een blad van de boom moet zitten. De functie tune() is hier nog een placeholder om de beste waarden voor deze parameters - die we later bepalen - daar in te stellen. We gebruiken 1.000 bomen c.q. versies van het model.
## Bouw het model: random forestrf_mod <-rand_forest(mtry =tune(), min_n =tune(), trees =1000) |>set_engine("ranger", num.threads = cores) |>set_mode("classification")
4.3 Maak de recipe
We maken een recipe voor het random forest model. We verwijderen de student ID en het collegejaar uit de data, omdat deze niet moet worden gebruikt in het model. Overige stappen zijn bij een random forest minder relevant in tegenstelling tot een regressiemodel.
## Maak de recipe: random forestrf_recipe <-recipe(Uitval ~ ., data = dfUitval_train) |>step_rm(ID, Collegejaar) ## Verwijder ID en Collegejaar uit het model## Toon de recipetidy(rf_recipe) |> knitr::kable()
number
operation
type
trained
skip
id
1
step
rm
FALSE
FALSE
rm_u6aow
4.4 Maak de workflow
We voegen het model en de recipe toe aan de workflow voor dit model.
## Maak de workflow: random forestrf_workflow <-workflow() |>add_model(rf_mod) |>add_recipe(rf_recipe)## Toon de workflowrf_workflow
══ Workflow ════════════════════════════════════════════════════════════════════
Preprocessor: Recipe
Model: rand_forest()
── Preprocessor ────────────────────────────────────────────────────────────────
1 Recipe Step
• step_rm()
── Model ───────────────────────────────────────────────────────────────────────
Random Forest Model Specification (classification)
Main Arguments:
mtry = tune()
trees = 1000
min_n = tune()
Engine-Specific Arguments:
num.threads = cores
Computational engine: ranger
4.5 Tune en train het model
We trainen en tunen het model in de workflow. We maken een grid met verschillende waarden voor de parameters mtry en min_n. We gebruiken de Area under the ROC Curve (AUC/ROC) als performance metric. Met de resultaten van de tuning kiezen we het beste model.
## Toon de parameters die getuned kunnen wordenrf_mod
Random Forest Model Specification (classification)
Main Arguments:
mtry = tune()
trees = 1000
min_n = tune()
Engine-Specific Arguments:
num.threads = cores
Computational engine: ranger
## Extraheer de parameters die getuned wordenextract_parameter_set_dials(rf_mod)
Collection of 2 parameters for tuning
identifier type object
mtry mtry nparam[?]
min_n min_n nparam[+]
Model parameters needing finalization:
# Randomly Selected Predictors ('mtry')
See `?dials::finalize` or `?dials::update.parameters` for more information.
## Bepaal de seedset.seed(2904)## Bouw het grid: random forestrf_res <- rf_workflow |>tune_grid(dfUitval_validation,grid =25,control =control_grid(save_pred =TRUE),metrics =metric_set(roc_auc))
i Creating pre-processing data to finalize unknown parameter: mtry
4.6 Kies het beste model
We evalueren de beste modellen en maken een ROC curve om de performance van het model te visualiseren. Vervolgens vergelijken we de prestaties van de modellen en kiezen we het beste model.
## Toon de beste modellenrf_res |>show_best(metric ="roc_auc", n =15) |>mutate(mean =round(mean, 6)) |> knitr::kable()
mtry
min_n
.metric
.estimator
mean
n
std_err
.config
9
12
roc_auc
binary
0.655394
1
NA
Preprocessor1_Model05
18
40
roc_auc
binary
0.653464
1
NA
Preprocessor1_Model16
11
29
roc_auc
binary
0.651693
1
NA
Preprocessor1_Model02
2
6
roc_auc
binary
0.651534
1
NA
Preprocessor1_Model22
5
36
roc_auc
binary
0.650870
1
NA
Preprocessor1_Model09
16
34
roc_auc
binary
0.650617
1
NA
Preprocessor1_Model07
4
26
roc_auc
binary
0.650585
1
NA
Preprocessor1_Model23
13
25
roc_auc
binary
0.650554
1
NA
Preprocessor1_Model25
3
8
roc_auc
binary
0.650364
1
NA
Preprocessor1_Model08
7
35
roc_auc
binary
0.650142
1
NA
Preprocessor1_Model14
20
10
roc_auc
binary
0.650142
1
NA
Preprocessor1_Model18
11
21
roc_auc
binary
0.649953
1
NA
Preprocessor1_Model15
15
20
roc_auc
binary
0.649541
1
NA
Preprocessor1_Model20
22
38
roc_auc
binary
0.649415
1
NA
Preprocessor1_Model01
17
24
roc_auc
binary
0.649320
1
NA
Preprocessor1_Model10
## Plot de resultatenautoplot(rf_res)
## Selecteer het beste modelrf_best <- rf_res |>select_best(metric ="roc_auc")rf_best|> knitr::kable()
mtry
min_n
.config
9
12
Preprocessor1_Model05
Toon code
## Verzamel de predictiesrf_res |>collect_predictions() |>head(10) |> knitr::kable()
## Bepaal de AUC van het beste modelrf_auc_highest <- rf_res |>collect_predictions(parameters = rf_best) |>roc_auc(Uitval, .pred_FALSE)## Voeg de naam van het model en de AUC toe dfModel_resultsdfModel_results <- dfModel_results |>add_row(model ="Random Forest", auc = rf_auc_highest$.estimate)
5 De uiteindelijke fit
In de laatste stap van deze analyse maken we het model definitief.
We testen het model op de testset en evalueren het model met metrieken en de Variable Importance Factor (VIF).
5.1 Combineer de AUC/ROC curves en kies het beste model
Eerst combineren we de AUC/ROC curves van de modellen om ze te vergelijken. We kiezen het beste model op basis van de hoogste AUC/ROC.
Toon code
## Combineer de AUC/ROC curves om de modellen te vergelijkenbind_rows(lr_auc, rf_auc) |>ggplot(aes(x =1- specificity, y = sensitivity, col = model)) +geom_path(lwd =1.5, alpha =0.8) +geom_abline(lty =3) +coord_equal() +scale_color_viridis_d(option ="plasma", end = .6)
Toon code
## Bepaal welke van de modellen het beste is op basis van de hoogste AUC/ROCdfModel_results <- dfModel_results |>mutate(number =row_number()) |>mutate(best =ifelse(auc ==max(auc), TRUE, FALSE)) |>arrange(number)## Bepaal het beste modelsBest_model <- dfModel_results$model[dfModel_results$best ==TRUE]sBest_model_auc <-round(dfModel_results$auc[dfModel_results$best ==TRUE], 4)
Het beste model is het Logistic Regression model met een AUC/ROC van 0.6567. Het Logistic Regression model heeft een AUC van 0.6567. Het Random Forest model heeft een AUC van 0.6554. We ronden de analyse verder af met het Logistic Regression model.
5.2 Maak het finale model
We maken het finale model op basis van de beste parameters die we hebben gevonden. Door in de engine bij importance de impurity op te geven, wordt het beste random forest model gekozen om de data definitief mee te classificeren.
## Test het ontwikkelde model op de testset## Bepaal de optimale parameters## Bouw het laatste modelif (sBest_model =="Logistic Regression") { last_lr_mod <-logistic_reg(penalty = lr_best$penalty,mixture =1) |>set_engine("glmnet") |>set_mode("classification")} elseif (sBest_model =="Random Forest") { last_rf_mod <-rand_forest(mtry = rf_best$mtry,min_n = rf_best$min_n,trees =1000) |>set_engine("ranger", num.threads = cores, importance ="impurity") |>set_mode("classification")}
5.3 Maak de workflow
We voegen het model toe aan de workflow en updaten de workflow met het finale model.
## Bouw de laatste workflow op basis van het laatste modelif (sBest_model =="Logistic Regression") { last_lr_workflow <- lr_workflow |>update_model(last_lr_mod)} elseif (sBest_model =="Random Forest") { last_rf_workflow <- rf_workflow |>update_model(last_rf_mod)}
5.4 Fit het finale model
We voeren de finale fit uit. De functie last_fit past het model toe op de validatieset.
## Voer de laatste fit uitset.seed(2904)if(sBest_model =="Logistic Regression") { last_fit <- last_lr_workflow |>last_fit(splits)} elseif(sBest_model =="Random Forest") { last_fit <- last_rf_workflow |>last_fit(splits)}
5.5 Evalueer het finale model: metrieken en vif
We evalueren het finale model op basis van 3 metrieken (accuraatheid, ROC/AUC en de Brier score = de Mean Squared Error) en de Variable Importance Factor (VIF). Uit de VIF is op te maken welke variabelen het meest bijdragen aan de voorspelling van de uitkomstvariabele.
## Verzamel de metriekenlast_fit |>collect_metrics() |>mutate(.estimate =round(.estimate, 4)) |> knitr::kable()
We weten nu welke variabelen van invloed zijn, maar niet hoe en in welke richting: dragen ze sterk bij of juist niet, verhogen of verlagen ze uitval? Om het model beter te begrijpen en te kunnen uitleggen, maken met behulp van het Dalex package een explainer. Dit package is ontwikkeld om beter uit te leggen welke variabelen van belang zijn en wat deze voor een effect hebben in een model. Een explainer is een model-onafhankelijke wrapper om het model heen en geeft inzicht in de voorspellingen van het model en de bijdrage van de variabelen aan de prognose. Een explainer maakt het verder mogelijk om modellen onderling te vergelijken en benchmarken.
6.1 Maak een explainer
Op basis van het tidymodels model extraheren we de informatie voor de explainer.
## Extraheer het fitted model en de workflowfitted_model <- last_fit |>extract_fit_parsnip()workflow <- last_fit |>extract_workflow()# Pas de Uitval variabele aan naar numeric (0/1), # zodat er een explainer van gemaakt kan wordendfOpleiding_inschrijvingen$Uitval <-as.numeric(dfOpleiding_inschrijvingen$Uitval) -1# Maak een explainerexplain_lm <- DALEX::explain(model = workflow,data = dfOpleiding_inschrijvingen,y = dfOpleiding_inschrijvingen$Uitval,label ="Linear Regression")
Preparation of a new explainer is initiated
-> model label : Linear Regression
-> data : 1818 rows 27 cols
-> target variable : 1818 values
-> predict function : yhat.workflow will be used ( default )
-> predicted values : No value for predict function target column. ( default )
-> model_info : package tidymodels , ver. 1.2.0 , task classification ( default )
-> predicted values : numerical, min = 0.09233298 , mean = 0.3963696 , max = 0.738768
-> residual function : difference between y and yhat ( default )
-> residuals : numerical, min = -0.7024894 , mean = 0.003520435 , max = 0.8469217
A new explainer has been created!
Toon code
library(iBreakDown)## Bereken de model parts op basis van de RMSEmp <-model_parts(explain_lm, loss_function = loss_root_mean_square)# Toon de resultatenplot(mp)
6.2 Pas het model toe op een willekeurige student
Nu we deze explainer hebben, kunnen we het model toepassen op een willekeurige student. We kiezen een student uit de dataset en passen het model toe op deze student. We kijken naar de voorspelling van het model en de bijdrage van de variabelen aan deze voorspelling. Dit geeft een verder inzicht in de werking van het model.
De intercept van het model (de basisvoorspelling) is 0,396: een kans op uitval van 39.6%. De bijdrage van de variabelen aan de voorspelling kan positief of negatief zijn. Een positieve bijdrage betekent dat de variabele de voorspelling verhoogt, een negatieve bijdrage betekent dat de voorspelling wordt verlaagd. De combinatie van variabelen bij deze student verlaagt de kans op uitval naar 0,254 (25,4%).
Toon code
source("01. Includes/Studentpersonas_opleidingen.R")## Loop over de persona'sfor (i in1:nrow(dfPersona)) { student_x <- dfPersona[i,] bd_lm <-predict_parts(explainer = explain_lm,new_observation = student_x,type ="break_down")# Convert to data frame for ggplot bd_lm_df <-as.data.frame(bd_lm) |>## Sorteer op basis van de positionarrange(desc(position)) |>## Maak het label aanmutate(label =paste0(as.character(round(contribution, 3)*100),"%")) |>mutate(label =case_when( variable %in%c("intercept","prediction") ~ label, sign ==1~paste0("+", label), sign ==-1~paste0(label), sign ==0~paste0("+", label),.default = label )) |>mutate(label =case_when( label =="0%"~"+0%",.default = label )) |>filter(sign !=0) |>mutate(start = dplyr::lag(cumulative, default =min(cumulative)),end = cumulative) ## Voeg een rij toe "Overige variabelen" bd_lm_df <- bd_lm_df |>add_row(variable ="+ Overige variabelen", contribution =0, variable_name =NA,variable_value =NA,cumulative = bd_lm_df$cumulative[bd_lm_df$position ==1],sign ="0", position =2,label ="+0%", start = bd_lm_df$cumulative[bd_lm_df$position ==1], end = bd_lm_df$cumulative[bd_lm_df$position ==1] ) |>arrange(position) |>mutate(position =row_number()) |>arrange(desc(position)) |>mutate(next_start =lead(start, default =NA),next_position =lead(position, default =NA)) |>mutate(start =case_when( variable =="intercept"~0, variable =="prediction"~ bd_lm_df$end[bd_lm_df$variable =="intercept"],.default = start)) |>mutate(sign =case_when( variable =="intercept"~"X",.default = sign)) |>mutate(label_color =case_when( sign =="1"~"#C8133B", sign =="-1"~"#466F9D",.default ="black")) |>rowwise() |>mutate(label_position =max(start, end)) |>ungroup() bd_lm_df$sign <-factor(bd_lm_df$sign, levels =c("1", "-1", "0", "X"), labels =c("Positief", "Negatief", "Geen", "X")) nUitval <-round(bd_lm_df$cumulative[bd_lm_df$variable =='prediction'], 3)*100 student_x_title <-paste0("Gemiddelde student | vooropleiding: ", student_x$Vooropleiding," | kans op uitval: ", nUitval, "%")# Bepaal de breaks voor de x-as y_breaks <-seq(0, 1, by =0.2) y_labels <-paste0(seq(0, 100, by =20), "%")# Watervalplot met ggplot p <-ggplot(bd_lm_df) +# Voeg horizontale lijnen toe om de 0.2geom_hline(yintercept = y_breaks, color ="#CBCBCB", linetype ="solid", linewidth =0.5) +geom_hline(yintercept =min(bd_lm_df$cumulative), color ="black", linetype ="dotted") +geom_rect(aes(xmin = position -0.4, xmax = position +0.4,ymin = start, ymax = end, fill = sign)) +geom_segment(aes(x = next_position -0.4, xend = position +0.4,y = end, yend = end), color ="darkgray") +coord_flip() +labs(title ="Breakdown van kans op Uitval",subtitle = student_x_title,x =NULL,y =NULL) +scale_fill_manual(values =c("Positief"="#C8133B", "Negatief"="#466F9D", "X"="#57606C")) +# Voeg tekstlabels toe voor de variabelengeom_text(aes(x = position, y = label_position, label = label),hjust =-0.1, size =4,color ="black") +# Pas het thema aan om de y-as labels weer te gevenscale_x_continuous(breaks = bd_lm_df$position, labels = bd_lm_df$variable) +scale_y_continuous(breaks = y_breaks, labels = y_labels, limits =c(0, 1)) +theme_minimal() +theme(legend.position ="none") +theme(axis.text.y =element_text(size =10),panel.grid.major =element_blank(), panel.grid.minor =element_blank() ) # Print de plotsuppressWarnings(print(p))}
6.3 Plot de ROC curve
Tot slot maken we een ROC curve om de prestaties van het definitieve model te visualiseren. De Sensitivity (True Positive Rate) en Specificity (True Negative Rate) worden hierin uitgezet. De Area under the ROC Curve (AUC/ROC) geeft de prestaties van het model weer. Het model scoort beter naarmate de AUC/ROC dichter bij de 1 ligt (in de linkerbovenhoek). Een AUC/ROC van 0,5 betekent dat het model niet beter presteert dan een willekeurige voorspelling.
## Toon de roc curvelast_fit |>collect_predictions() |>roc_curve(Uitval, .pred_FALSE) |>autoplot()
7 Conclusies
7.1 Het beste prognosemodel voor deze opleiding
Het beste prognosemodel blijkt het Logistic Regression model te zijn.
Van de prognosemodellen die we hebben ontwikkeld om uitval na 1 jaar te voorspellen, had het Logistic Regression model de hoogste AUC/ROC waarde (0.6567).
7.2 Mate van accuraatheid en lift
Een prognosemodel moet minimaal beter presteren dan een base-model om waarde toe te voegen. Het base-model neemt de grootste klasse van de gemiddelde uitval na 1 jaar van de afgelopen jaren als basis. Stel we zouden tegen alle studenten zeggen dat ze hun studie gaan halen, dan is de mate van accuratesse gelijk aan dit base-model. Dit base-model is dus altijd hoger dan de 50% lijn van de AUC/ROC curve, tenzij het base-model toevallig precies 50% is
De mate van accuraatheid van de toepassing van het model is vrij laag (65.48%).
Base-model: 60.01% – Voor deze opleiding bereken we het base-model als volgt. Van alle studenten viel 39.99% uit. De grootste klasse (en de accuratesse) van het base-model is daarmee (100% - 39.99% = ) 60.01% die niet uitviel.
Accuratesse prognose: 65.48% – Het model voorspelt Uitval na 1 jaar met een accuratesse van 65.48%.
Lift: 5.47% – Het model scoort in de huidige opbouw met een verschil van 5.47% (de lift) wat beter dan de accuraatheid van het base-model (60.01%).
7.3 Factoren
Vooralsnog geen nadere bijzonderheden.
Verantwoording
Deze analyse maakt deel uit van het onderzoek naar kansengelijkheid van het lectoraat Learning Technology & Analytics van De Haagse Hogeschool: No Fairness without Awareness | Het rapport is door het lectoraat ontwikkeld in Quarto 1.4.549. | Template versie: 0.9.1.9000